General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



NASA Contractor Report 156886 (I) 


(MASA'^E-‘156686'’?ol-1) JIOSS lIGOfilTMH 

CUfifiEiT H1PPII6, 
VOLuas 1 tAnalytic Sciences Cocp.) 45 n 

HC »03/HF A01 CSCl 08C 


fll/46 


M82-30623 


Oncias 

28727 


NOSS Algorithm Specifications for 
Ocean Current Mapping 


Volume I 


James V. White 



NASA 

National Aeronautics and 
Space Administration 

Goddard Space HigM Center 

Wallops Flight Center 
Wallops Island. Virginia 23337 


NASA Contractor Report 156886 (I) 


NOSS Algorithm Specifications for 
Ocean Current Mapping 

Volume I 


James V. White 

The Analytic Sciences Corporation 
One Jacob Way 

Reading. Massachusetts 01867 


Prepared Under Contract No. NAS6-3163 


NASA 

National Aeronautics and 
Space Administration 

Goddard Space RieM Canter 

Wallops Flight Center 
Wallops Island. Virginia 23337 



PRECEDING PAGE BLANK NOT FILMED 


TABLE OF CONTENTS 

Page 

No. 


List of Figures iv 


1. OVERVIEW 1-1 

1.1 Introduction 1-1 

1.2 Algorithm Descriptions 1-1 

1.3 Generic Ocean-Current Signatures 1*3 

2. AUTOREGRESSIVE MODELING ALGORITHM 2-1 

2.1 Theory 2-1 

2.2 Augmented COVAR Algorithm 2-2 

2.3 Specification for the Autoregressive 

Modeling Algorithm 2-5 

3. MATCHED-FILTER DESIGN ALGORITHM 3-1 

3.1 Theory ' 3-1 

3.2 Specifications for the Matched-Filter 

Design Algorithm 3-6 

A. MATCHED-FILTER CONVOLUTION ALGORITHM A-1 

A.l Theory A-1 

A. 2 Specifications for the Matched-Filter 

Convolution Algorithm A-1 

5. THRESHOLD DETECTION ALGORITHM 5-1 

5.1 Theory 5-1 

5.2 Specifications for the Threshold 

Detection Algorithm 5-3 

6. GEOSTROPHIC-VELOCITY ESTIMATION ALGORITHM 6-1 

6.1 Theory 6-1 

6.2 Specifications for the Geostrophic-Velocity 

Estimation Algorithm 6-3 

7. SUMMARY 7-1 

REFERENCES R-1 


iii 


LIST OF FIGURES 


Figure 


Page 

No. 


No. 

1.2-1 

Structure of Dete-Adaptive Current 
Detection Algorithm 

1-2 

1.3-i 

Gaussian Cold-Ring Signature, Depth = 0.5 m, 
Width = 150 km 

1-5 

1.3-2 

Tangential Current-Velocity Distributions in 
Gaussian Cold Ring 

1-6 

1.3-3 

Geometry of Geostrophic Current and 
Satellite Subtrack 

1-7 

1.3-A 

Dynamic Sea-Surface Height Signature 

1-9 

1.3-5 

Geostrophic Velocity Profile 

1-9 

3.1-1 

Matched-Filter Detector 

3-2 


iv 


1 . 


OVERVIEW 


1 . 1 INTRODUCTION 

This chapter gives an overview of the algorithms for 
detecting ocean currents and estimating geostrophic velocities 
with satellite altimeter data. Chapters 2 through 6 describe 
the theory of operation and the specifications for each algorithm. 


1 . 2 ALGORITHM DESCRIPTIONS 


The current-detection and velocity-estimation algo- 
rithms process single tracks of residual satellite altimeter 
data and yield the following outputs: 


• Detected locations of specified ocean- 
current signatures along the satellite 
subtrack 

• Estimated amplitudes of the detected 
signatures 

• Estimated rms errors for the locations 
and amplitudes of detected signatures 

• Estimated cross-track component of the 
boundary -current geostrophic velocity 
and an rms error bound for the estimate 

• Expected number of false alarms. 


The residual altimeter data are inputs to the algorithms 
and are computed from raw altimeter data in three steps by 
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• Applying corrections for known error 
sources 

e Interpolating the data through intervals 
in which the data are in serious error 
(e.g., outliers) 

a / ibtracting an estiaated graviaetric 

geoid profile along the satellite subtrack. 

The resulting residual data are noisy measurenents of 
the dynamic sea-surface height. The characteristics of the 
noise in these depend on the noise in the raw altimeter data 
and on the accuracies of the error corrections and the geoid 
profiles. The detection algorithm exploits both the statistical 
properties of the noise in the residual data and the known aver- 
age properties of ocean-current signatures in the altimeter data. 
For specified models of the noise and oceanographic signature* 
the algorithm maximizes the probability of detection at a spe- 
cified probability of false alarm and minimizes the rms errors 
in the estimated current signature parameters. 


As depicted in Fig. 1.2-1, the detection algorithm 
consists of four subalgorithms that perform separate functions. 



Figure 1.2-1 Structure of Data-Adaptive Current 

Detection Algorithm 
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• AUTOREGRESSIVE MODELING - The residuel 

date are anelyzed to deternine a stochas- 
tic autoregressive (AR) model for the 
process that generated the data. The 
order of the AR model is selected to 
minimize the Akaike information criterion. 

e MATCHp-FILTER DESIGN - The AR model, 
together with a user-specified generic 
ocean-current signature, are used to 
compute the impulse response of the op- 
timal matched filter lor detecting and 
locating the generic signature in the 
noisy residual data. 

e MATQIED FILTER - The impulse response of 
the matched filter is convolved with the 
residual altimeter data to compute a 
sequence of sufficient statistics for 
the threshold detector. 

e THRESHO^ DETECTOR - The detector compares 
the. sufficient statistics with a threshold 
value that is chosen to yield a specified 
false-alarm rate. A detection occurs 
when the statistic exceeds the threshold. 
The estimated location of the detected 
signature is given by the location of 
the local maximum of the statistic. 


1.3 GENERIC OCEAN-CURRENT SIGNATURES 

This section describes two parametric families of 
ocean-current signatures. The first family is used for design- 
ing matched filters to detect warm-core and cold-core current 
rings. The second family is intended for detecting boundary 
currents, such as the Gulf Stream, and for estimating geo- 
strophic current velocities. 

Ring-Current Signatures - A family of generic alti- 
metric signatures is described for modeling the dynamic sea- 
surface features caused by cold-core and warm-core current 
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rings. The ses-surfece height H(x) at radial position x with 
respect to the ring's center is Modeled as 

H(x) s - D exp <-9.21 <x/W)^) (1.3-1) 

D = signature depth 
V = signature width 

D is positive for cold rings and negative for warn rings. This 
parametric model has a simple mathematical form and appears to 
be in reasonable agreement with available data on ring signa- 
tures <e.g., Refs. 1-3). 

The width V is defined as the diameter at which the 
signature is 10 percent of its central value: 

H(W/2) = H(0)/10 (1.3-2) 

Equation 1.3-1 has the form of a Gaussian probability 
density. Therefore, these are referred to as Gaussian ring 
signatures. An example of a Gaussian ring signatun* is shown 
in Fig. 1.3-1, where the central depth is 0.5 meter and the 
width is 150 km. 

The tangential current -velc/city distribution implied 
by a ring signature may be computed by setting the radial slope 
of the sea surface equal to the sum of the horizontal Coriolis 
acceleration and the centrifugal acceleration divided by the 
acceleration of gravity: 

dH(x) . f v(x) ♦ v^(x)/x 
dx " g 


(1.3-3) 



QA!J88IAN COLO-HINO MOOfL m-vmis. 



Figure r.3*l Gaussian CoId*Ring Signature, 

Depth = 0.5 B, Width = 150 km 

£ = 20 sin# = Coriolis parameter 
0 s earth's rotational velocity 
# = latitude 

v(x) = tangential current velocity 
g - acceleration of gravity. 

The geostrophic velocity component is 

(1.3-4) 


Solving Eq. 1.3-3 for the total current velocity v(x) yields 

(1.3-5) 


v(x) = 


J 

Kl + ~A] r - 1 


ITT 
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For « Gaussian 0.5-a 150-loi ring signature at 45-degrees 
latitude. Eqs. 1.3-4 and 1.3-5 yield the velocity distributions 
shown in Fig. 1.3-2. The geos trophic approxiaation is seen to 
over-estimate the maximum velocity by approximately 0.1 m/s <14%). 


TANGENTIAL CURRENT VELOCITY 



Figure 1.3-2 Tangential Current-Velocity Distributions in 

Gaussian Cold Ring. Crosses s geos trophic 
approximation; circles = geostrophic approxi- 
mation with centrifugal correction. 


Boundary-Current Signatures - A family of generic al- 
timetric signatures for boundary currents is defined with the 
aid of Fig. 1.3-3, which depicts a satellite subtrack crossing 
a current at angle 6. At position x along the subtrack, the 
dynamic height H(x) is modeled with the hyperbolic tangent 
function. 

H(x) = -(V2) tanh(3 x sin0/V^) (1.3-6) 
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Figure 1.3-3 Geonetry of Geostrophic Current and 

Satellite Subtrack 

A = amplitude of dynamic height change 
6 = track angle with respect to current velocity 
V = width of current (90% height change) 

The along-track slope of the signature is 
dH(x) _ 3 A sin0 »_2r3 x sinB 

- - 2 w, ■ [—r — 

C L c 

For the coordinate system in Fig. 1.3-3» this slope is related 
to the cross- track component V (x) of the geostrophic velocity 
as follows 


(1.3-7) 




For tracks chat cut across the current* Che geostrophic 
velocity profile V^(x) along the subtrack is proportional to 
Che cross-track velocity V^(x) 


Vg(x) * V^(x)/sine 

For the signature slope given by Eq. 1.3-7* 
strophic velocity profile is therefore 



(1.3-9) 
the geo- 


( 1 . 3 - 10 ) 


The signature aaplitude parameter A is proportional 
to the aaxinusi geostrophic velocity V^(0) 

2 f W 

A = ■ V.(0) (1.3-11) 

3 g g 


For tracks that intersect the current, the signature width 
is proportional to the current width 

W. * W Veins (1.3-12) 

s c 

Typical »odel parasieters for the Gulf Stress in the 
western North Atlantic are an aaplitude of A « 1 a saxisnui 
geostrophic velocity of V^(0) = 2 s/s* and a latitude of 
S > 45 deg. Froa Eq. 1.3-11 the current's width is > 71 ka. 
For a ncMiinal track crossing angle of 60 deg* the along- track 
width of the signature is * 82 ka. Figure 1.3-4 shows the 
dynaaic sea-surface height signature for these paraaeter values* 
while Fig. 1.3-5 depicts the geostrophic velocity profile ia- 
plied by the height signature. 
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Figure 1.3-5 Geostrophic Velocity Profile 


1-9 




2 . 


AUTOREGRESSIVE MODELING ALGORITHM 


2.1 THEORY 


An autoregressive (AR) nodel of order p for a tise 

series {x(t>; t 1»2 n} is a difference equation driven 

by white noise w^(t) 


x(t) = Cj x(t-l) + x(t-2) ♦ ... + Cp x(t-p) ♦ Wp(t) 

t = p+1, p+2, . . . , n 


( 2 . 1 - 1 ) 


To identify the best AR nodel for the underlying random 
process that generated the x(t) data, a family of AR models is 
considered. Each member of the family corresponds to a different 
model order 


p = 0,1,2, .. . ,pmax 


When modeling residual altimeter data, sampled at a 1-Hz rate, 
a reasonable choice for the maximum order is 

pmax = integer (n/20) 

For each order (p = 0,1,2 pmax) in turn, the AR 

coefficients ,€ 2 , . • • ,Cp (C^ = 1) are selected to minimize the 
sample noise variance 


a2(p) = _1 V w^(ii) 

n-pmax ,=ptaxn P 


( 2 . 1 - 2 ) 
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This variance is then used to compute the Akaike infomation 
criterion (Refs. 5-7) 

AIC(p) = n In(o^) ♦ 2p (2.X-3) 


The particular order p that ninimites AIC is identified; the 
corresponding AR paraneters 


then define the particular AR nodel (for the x(t) process) 
that is best supported by the available data. 


There are many known algorithms for c<Miputing the AR 
coefficients to minimize the sample variance . in Eq. 2.1-2. An 
effective algorithm for the present application is based on an 
augmented version of the COVAR algorithm (Ref. 4). which is 
described in the following section. 


2.2 AUGMERTED COVAR ALGORITHM 


This section describes an augmented version of the 
COVAR algorithm (Ref. 4), which solves the following problem 
by Cholesky factorization. 


GIVEN: 1. Data sequence {x(0),x(l) 


2. Integer M 


,x(N-D) 


FIND: 


1. Coefficients {a. ,a 2 > . . . »Sw) that minimize 
the suffl-squared‘'’AR^re8iduals 


or = 2 + E x(n-k)l^ (2,2-1) 

n=M k=l 

2. Minimized value of or 
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SOLUTION; 


I. Define 


N-1 


'lit " x(n-i) x(n-k); k = 1,2 M 

1 s 1,2,. ...M 


( 2 . 2 - 2 ) 


2. Solve the following M equations for 

" • «ok*' ^ ’ 1,2,...,M (2.2-3) 


In the following description of the augnented COVAR algorithm, 
asterisks are used to indicate additions to the original algo- 
rithm in Ref. 4. 

INPUTS: N, (x(t); t=0,l, . . . ,N-1} , M 

N > number of time-series data 
x(t) • datum at time t 

M = maximum autoregressive order 

OUTPUTS: M. (C(k); k = 0,1,2 M) , {HS(k); k = 0,1,2,...,M} 

M = maximum AR order 
C(k) = k^^ AR coefficient 

* MS(k) = mean- square value of AR residuals for model 
of order k 


ALGORITHM: 

1. Compute Cqq, c^^q, and using Eq. 2.2-2 

2. Initialize the following parameters 

®00 “ ^ (2.2-4) 

(2.2-5) 


®0 " ^00 
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3. 


* MS(0) = «q/(N-M) 

*10 “ ^ 

*11 ‘ ^1 
^01 = 1 
^0 ‘ ®11 

“1 ' “0 - '‘l^O 

* MS(1) = a^/(N-M) 

Recursively stepping m (for m = 1,2 
'■+ 1.0 2 - 2*2 
'■+l,k = '..k-l * x(M-k) 

- x(N-a-l) x(N-k); k = 1,2 

, n+1 

^■n ' f; '»+l.j '’nj* ” = 0-2 



^m,iB+l 




'ai ”ij» 


j = 1,2.. 


1 


B ^l 


( 2 . 2 - 6 ) 

(2.2-7) 

( 2 . 2 - 8 ) 

(2.2-9) 

( 2 . 2 - 10 ) 

( 2 . 2 - 11 ) 

( 2 . 2 - 12 ) 

(2.2-13) 

. ..M-1), coapute: 


...a+l (2.2-14) 

...a-1 (2.2-15) 

.a (2.2-16) 
(2.2-17) 
(2.2-18) 


‘‘■+1 ' * 'nfH,! *ni 

*■+1,0 * ^ 

*■+1,1 * *mi * *^0+1 **mi' ^ 

*m"*'l ,a+l ” ^m+1 

"■+1 ' “m * ‘‘■+1 ^ 


(2.2-19) 


( 2 . 2 - 20 ) 
,a (2.2-21) 
( 2 . 2 - 22 ) 
(2.2-23) 
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* MS(m+l) = 

STEP m OF RECURSION IS COMPLETED 
4. Termination (end of step M-l) 

^ ~ »M 

a - a 

m 

* C(0) = 1 

* C(k) = -aj^; k = 1»2, . . . ,M 


(2.2-24) 


(2.2-25) 

(2.2-26) 

(2.2-27) 

(2.2-28) 


The augmented COVAR algorithm is called as a subroutine 
in the AR modeling algorithm (ACOVAR) used for ocean-current 
detection, which is specified below. 


2.3 SPECIFICATION FOR THE AUTOREGRESSIVE MODELING ALGORITHM 


For ocean-current detection, the ACOVAR algorithm is 
used for autoregressive modeling. ACOVAR uses the augmented 
COVAR algorithm (Section 2.2) as a subroutine. Formal speci- 
fications for ACOVAR are given in the following. 


NAME: ACOVAR 

PURPOSE: Compute the parameters of an optimal autoregressive 
(AR) model for one track of residual altimeter data. 

INPUTS: N, (D(k); k = 1,2 N} 

N = number of residual altimeter data 

D(k) = k^^ sample in the time series of residual altimeter 
data along one satellite subtrack 


OUTPUTS: P, {C(k); k = 0,1,...,P}, VAR 
P = order of selected AR model 
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C(k) = AR-nodel coefficient 
VAR = nean- square value of AR residuals 


ALGORITHM: 


1. Use the augmented COVAR algorithm to compute the AR 
coefficients C<k) and the mean-square residuals MS(k) 
(for k = 0,1,2....,M) by using the following assignments 

x(k) = D(k+1); k = 0.1 N-1 (2.3-1) 

M = INTEGER (N/20) (2.3-2) 

Set PMAX = M. (2.3-3) 


2 . 


Compute the Akaike information criterion {AlC(k) ; 
k = 0 ,1 . . . . .PMAX} for each order of AR model 

AlC(k) = N ln(MS(k)) + 2k 


k = 0,1,. . . ,PMAX 


(2.3-4) 


3. Determine the smallest J such that 


0 < J < PMAX 

• • (2.3-5) 

AIC(J) i AIC(I); I = 0.1,..., PMAX 

4. If (J = PMAX) or (J = 0) then 

P = J (2.3-6) 

VAR = MS(J) (2.3-7) 

GOTO STEP 5 
If (J < PMAX) then 

GOTO STEP 1 BUT USE M=J (2.3-8) 

5. OUTPUT P, (C(k); k = 0,1 P} , VAR 

6. END 
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MATCHED-FILTER DESIGN ALGORITHM 


3.1 THEORY 

The theory of optimal matched filters for detecting 
deterministic signatures in additive colored noise is discussed 
in several text books (e.g., Refs. 8, 9, 11). The key results 
of that theory for ocean-current detection are summarized in 
the following. 

The problem of detecting ocean-current signatures in 
residual altimeter data is formalized as follows. 

GIVEN: D(t) = time series of residual altimeter data 

m(t) = ocean-current signature time series 

N(t) = stationary Gaussian noise model for residual 
altimeter data that are free of m(t) 

T = specified time (location) in the data D(t) 

Ag = unknown signature amplitude scale factor 

Rp = hypothesis that D(t) = N(t) + A_m(t-T) with 

Hq = hypothesis that D(t) = N(t) 

FIND: An optimal decision rule for correctly choosing between 

hypotheses Hq and H..; and an optimal estimate of the 
amplitude A^^when Rj. is chosen. 

OPTIMALITY: Maximize the probability of correct detection 

for a specified probability of false alarm. 

SOLUTION: Compute the likelihood ratio 
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Likelihood of D(t) under Hj, 

" Likelihood of D(t) under Hq 

Select Hp when LR > threshold value. 

Select Hq when LR 1 threshold value. 

As depicted in Fig. 3.1-lt the optimal decision rule can be 
efficiently implemented by processing the residual altimeter 
data D(t) with one matched filter and a threshold detector to 
test Hp against Hq for all possible values of T. Once a de- 
tection is made (i.e., Hp is selected), the best estimates of 
the location T and the amplitude scale factor are easily 
computed from the matched- filter output. 


nmm 


RECOMMENDED 8CAUNG Of FILTER OUPUT 
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Figure 3.1-1 Matched-Filter Detector 

The optimal matched filter for long data sets is a 
convolution operator having the frequency response H(F) and 
the impulse response h<t) (e.g.. Ref. 11, pp. 325-329) 
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(3.1-1) 

( 3 . 1 - 2 ) 


M(F) = Fourier transform of the ocean-current signature- 
m(t) 

Sjjjj(F) = power spectrum of the res iual altimetry N(t) 

The Fourier transform of m(t) is defined as 

M(F) = S m(t) 

t=-« 

M(F) dF 

When the residual altimeter noise model N(t) is auto- 
regressive (AR) , the optimal matched filter can be implemented 
as a finite- impulse- response (FIR) filter. This means that 
the matched- filter impulse response h(t) has finite support, 
i.e., h(t) = 0 for t < and t > for finite and 

^max' 



(3.1-3) 

(3.1-4) 
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h(t) H(F) dF 


H(F) = 

F = frequency (cycles/sample) 


The AR model for the noise N(t) is a difference equa- 
tion of order p driven by white noise W(t) 

N(t) = N(t-l) + C 2 N(t-2) + ... + Cp N(t-p) + W(t) (3.1-5) 
t — ...—1,0,1... 

Mean(W(t>) = 0; Variance(W(t)) = 

From linear-system theory (Refs. 10 and 11), the power 
spectnun Sj^j|(F) of N(t) is 
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* a(F)°c ; -F) (3.1-6) 

G(F) * 1 - i C,J ,-12iiFk (3.1-7) 


From Eqs. 3.1*2 and 3.1-6, the matched- filter frequency response 
may be expressed as 



o-2g(F) G(-F) M(-F) 


(3.1-8) 


The Fourier transform of Eq. 3.1-8 yields the following 
expression for the impulse response h(t) of the matched filter 


h(t) = o'^(t) * g(-t) * «(-t) (3.1-9) 
G(F) e^^nFt (3.1-10) 

g(t) s 0; t < 0 (3.1-11) 
g(0) = 1 (3.1-12) 
g(l) = -Cj (3.1-13) 


g(p) = -Cp (3.1-14) 

g(t) = 0; t > p (3.1-15) 

Since the convolutions in Eq. 3.1-9 contain only a finite number 
of non-zero terms when m(t) has finite support, the impulse 
response h(t) also has finite support. 

The rms signal -to-noise ratio achieved by the matched 
filter is defined as 
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« Peak Filter Output Due to Sl£nature m(t) 
~ Rms Filter Output Due to Noise Vi(t) 


(3.1-16) 


The SNR of an optimal matched filter can be computed 
with the formula 


SNR = J2- h(j) m(-j> (3.1-17) 

lj = -* 

where only a finite number of terms contribute. The rms value 
of the noise in the filter output is numerically equal to SNR 
when the filter is optimized 

Rms Noise in Filter Output = SNR (3.1-18) 

The peak filter output value due to the signature m(t) is also 
expressible in teras of SNR: 

2 

Peak Filter Output Due to Signature m(t) = SNR (3.1-19) 

Since the SNR equals the rms value of the modeled 
noise in the filtered output, it is convenient to scale the 
filter output by the factor 1/SNR as shown in Fig. 3.1-1. 

This yields a test statistic Y(t) that contains a random com- 
ponent having a standard deviation of unity. 

The mean output frequency F^^^ of the filter is a n\im- 
ber that measures the average rate at which the filter output 
changes sign 

F = Half the Average Rate of Zero Crossings (3.1-20) 
® of Noise in Filter Output 
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By using results in Ref. 10, p. 199, it My be shown for the 
Geussien noise in the filter output that 

^ cos“^(x/SNR^l (cycles/ssBple) (3.1-21) 

m 

X = hik) B<l-k) 

k»-» 


Equation 3.1-21 is derived in the Appendix of Ref. 12. 


3.2 SPECIFICATIONS FOR THE HATCHED- FILTER DESI(»< ALGORITHM 


For a specified AR noise Bodel and a specified ocean- 
current signature, the DESIGNMF algorithm is used to design 
the optimal matched filter. Formal specifications for this 
algorithm are given in the following. 


NAME: DESIGNMF 


PURPOSE: Compute the impulse response of the matched filter, 
the signal- to-noise ratio, and the mean output 
frequency . 


INPUTS: P, {C(k); k = 0,1,...,PJ, VAR, NM NH, 

{M(k); k = -NH,-NM’^1,...,NM}, SNR, FM 

P = order of the AR model for the residual altimeter 

data 

C(k) s k^^ coefficient of the AR model 

VAR = mean-square value of the AR* residuals 

NM = half -width of the time series containing the 

oceanographic signature to be detected (support 
of signature contains 2NM 1 points) 

M(k) = k^^ sample in the oceanographic signature time 
series 
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'Oi = half-width of the Batched- filter iapulae response 

(support of impulse response contains 2NH * 1 points) 

SNR = rms signal-to-noise ratio 

FM = mean output frequency (cycles/saaple) 


OUTPUTS: NH. {H(k); k * -NH.-NH+l NH} 

NH = half-width of the Batched- filter impulse 
response 

H(k) = k^^ sample in the matched- filter impulse 
response 


ALGORITHM: 

1. Compute A(k), the convolution of G(k) and G(-k): 

G(k) = 0, k < 0 (3.2-1) 

G(a) = 1; k = 0 (3.2-2) 

G(k) = -C(k); k * 1.2 P (3.2-3) 

G(k) s 0; k > P (3.2-4) 

A(k) = 2 C(j) G<J-k)i k = -P.-P+l P 

(3.2-5) 

2. Compute B(k). the convolution of A(k) and H'(-k): 

M’(k) = 0; k < -NM (3.2-6) 

M'(k) = M(k); k = -NM.-NM+l , . . . .NM (3.2-7) 
M’(k) = 0; k > NM (3.2-8) 

B(k) = I; A(j) M'(j-k); k = -NH.-NH+l , . . . ,NH 

(3.2-9) 

3. Compute H(k), the impulse response of the matched 
filter, by scaling B(k): 

H(k) =VAR'%(k); k = -NH.-NH+1 NH (3.2-10) 
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4. CoMpute SNR» the rmt eignel-to-noite retio: 

NM 

SNR * ( 2 (3.2-11) 

ka>-NN 

5. Conpute FH, the seen output frequency: 

FH s (2PI)*^ cos’^(x/SNR^) (cyclee/ea»ple) 

(3.2-12) 

NM 

X = ^ H(k) M(l-k) (3.2-13) 

k=-NM 

PI = 3.14159... 

6. Output: ‘nH, (H(k); k = -NH.-NH+l . . . . ,NH) 

7. End 
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MATCHED-FILTER CONVOLUTIOM ALGC«ITHM 


A. 


4.1 THEORY 

The optimal matched filter, depicted in Fig. 3.1-1, 
processes the residual altimeter data D(t) to produce a test 
statistic Y(t) for the threshold detector. The scaled output 
of the matched filter is computed through the following convo- 
lution of the data D(t) with the filter's impulse response h(t) 


Y(t) = 5515 E hOt) D(t-k) 


(A. 1-1) 


Since both h(t) and D(t) are finite sequences, the >convolution 
in Eq. 4.1-1 has only a finite number of terms. 


4.2 SPECIFICATIONS FOR THE MATCHED- FILTER CONVOLUTION ALGORITHM 

The matched filter is implemented as a convolution in 
algorithm CONVOLVE. The formal specifications for this algorithm 
are given in the following. 

NAME: CONVOLVE 

PURPOSE: Convolve the residual altimeter data with the matched- 
filter impulse response and scale the output. 

INPUTS: NH, {H(k) ; k = -NH,-NH+1 NH} , N, 

{D(k); k = 1,2,. . . ,N} , SNR 

NH = half-width of the matched-filter impulse response 


4-1 


H(k) = sample of the aetched- filter iapulse response 
N = musber of data in the residual altimetry tiae series 
D(k) =■ k^^ sasiple In the residual altimetry time series 
SNR = rms signal* to>noise ratio 

OUTPUTS: N, {Y(k); k = 1.2 N} 

N =: number of samples in the scaled output of the 
matched filter 

Y(k) = k^^ sample in the scaled output of the matched filter 
ALGORITHM: 


1. 

Ccmipute Y(k). the scaled convolution of H(k) 

and L' (k) 


D'(k) 

= 0; k < 1 

(4.2-1) 


D'(k) 

= D(k); k = 1.2.. ...N 

(4.2-2) 


D*(k) 

S 0; k > N 

(4.2-3) 



, NH 



Y(k) 

s SNR*-^ £ H(j) D’(k-j); 

(4.2-4) 



j=-NH 




k = NH+l.NH+2.. . . ,N-NH 



Y(k) 

= 0; k = 1,2,. . . ,NH 

(4.2-5) 


Y(k) 

= 0; k = N-NH+l.N-NH+2,. . . ,N 

(4.2-6) 

2. 

Output: N 

. {Y(k>; k = 1.2... ,N} 


3. 

End 
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5. 


THRESHOLD DETECTION ALGORITHM 


5.1 THEORY 

The threshold detector compares the sequence of scaled 
test statistics Y<t) from the matched filter against a detection 
threshold. When Y(t) exceeds the threshold, an alarm is said 
to occur. These alarms are classified into three categories: 

• Correct Detection caused by occurrences 
of modeled cur rint signatures in the 
residual altimeter data 

• False Alarms caused by random excursions 
of -the modeled noise in the residual 
altimeter data 

• Unmodeled Detections caused by unmodeled 
^current signatures or uMtodeled noise In 
the residual altimeter data. 

The statistics of correct detections and false alarms 
are computed for the specific ocean-current signature and the 
specific noise model for which the filter was designed. On 
the basis of these statistics, the expected average performance 
of the detector is predicted, and the detection threshold is 
adjusted for a desired tradeoff between detection probabilities 
and false-alarm rates. 

Formulas are listed below for computing the follow- 
ing performance statistics of the detector: the probability 

of false alarms; the average false-alarm rate, the maximum- 
likelihood estimates of signature location and amplitude (and 
their rms accuracies), and the probability of detecting a sig- 
nature with a prescribed amplitude. 
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Probability of False Alara - Let TH denote the detec- 
tion threshold. The probability that the noise component 
alone in Y(t) will exceed TH is given by the standardized nonal 
probability distribution function: 

P^ = Prob{Y(t) > TH; noise alone} (5.1-1) 

P^ = Q(TH) sjQ exp(-xV2) dx (5.1-2) 

The detection threshold that yields a prescribed probability 
of false alam is 


TH = Q'^(P£) (5.1-3) 

Averaxe False-Alam Rate - The threshold detector 
processes data fifoa individual tracks of altiaeter data» and 
it is often reasonable to set the detection threshold so that 
a specified nunber of false alaras is expected to occur per 
unit distance along the track (expressed in the units of alarms 
per data sample). This false-alarm rate (FAR) is computed as 

FAR = (F.) exp(-THV2) (alarms/sample) (5.1-4) 

m 

F • mean output frequency of modeled noise 
* in filter output 

Equation 5.1-4 is derived from results in Ref. 10, p. 492. The 
expected number of false alarms (EN) along a track of data 
having N samples of the test statistic Y(t) is 


EN s N-FAR (5,1 t5) 

Detection Threshold - The detection threshold TH is 
chosen to yield a specified false-alarm rate FAR 
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TH = /-i ln(FAR/F^) 


(5.1-6) 


Equation 5.1-6 is obtained by solving Eq. 5.1-4 for TH. 

Detected Signature Location and Aaplitude - The best 
estimate of the location of a detected signature is the value 
of t for which Y(t) achieves its local maximum value above the 
threshold TH. Let t^ denote this estimate of signature location. 
The Cramer-Rao (C-R) lower bound on the rms error of this esti- 
mate depends on the filter's maximum scaled output Y(t^): 

C-R Lower Bound = F \(t ) (5.1-7) 


The maximum- likelihood estimate of the signature amplitude 

scale factor A is 
s 



Y(t^) 

“5R!T 


(5.1-8) 


and the rms (one-sigma) error in this estimate is 1/SNR. The 

A 

best estimate of the detected signature is then Ag*m(t-t^). 

Probability of Detection - The probability of detect- 
ing the signature A m(t) is 

s 

Pd<As) = Q(TH - Ajj-SNR) (5.1*9) 

5.2 SPECIFICATIONS FOR THE THRESHOLD DETECTION ALGORITHM 

The threshold detector is implemented by algorithm 
DETECTOR. The formal specifications for this algorithm are 
given in the following. 
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NAME: 


DETECTOR 


PURPOSE: C<»p«re the scaled output of the Batched filter 

against a detection threshold; locate and count all 
alanss; compute a signature amplitude scale factor 
for each detection. 


INPUTS: N, {Y<k); k = 1,2 N} . NH. SNR. FM. FAR 

N = number of samples in scaled output of matched filter 

Y(k) = k^^ sample in scaled output of matched filter 

NH = half-width of matched- filter impulse response 

SNR = rms signal- to-noise ratio 

FH = mean output frequency (cycles/sample) 

FAR = average false-alarm rate (alarms/sample) 

(FAR must be < FM) 


OUTPUTS: NALARM, {A(k); k = 1,2 NALARMJ . 

{L(k); k = 1,2,. .. .NALARM} , EN, 

{CR(k); k = 1,2,. .. .NALARM} , AE 

NALARM = number of alarms 

A<k) = signature amplitude scale factor for k^^ 
detection 

L(k) = location of the k^^ detection (sample number) 
EN = expected number of false alarms 

CR(k) = Cramer-Rao lower bound on rms location 
error of k^^ detection (samples) 

AE = standard deviation of the errors in the 

estimates of the signature amplitude scale 
factors 


ALGORITHM: 


1. Compute detection threshold TH 


TH = V2 ln(PH/FAft) (5.2-1) 

2. Count number of elerms (NALARM) » record the location 

L(k) of the positive-going threshold crossing, end 
mark locations in array A(M) of all samples Y(k) that 


exceed the threshold 

NALARM = 0 (5.2-2) 

FLAG = 0 (5.2-3) 

FOR K = NH + 1 to N - NH 

IF FLAG = 0 AND Y(K) > TH THEN 

NALARM : NALARM >>> 1 (5.2-4) 

FLAG = 1 (5.2-5) 

L(NALARM) = K (5.2-6) 

IF Y(K) > TH THEN A(K) = 1 (5.2-7) 

IF Y(K) i TH THEN FLAG = 0 (5.2-8) 

NEXT K 


3. Skip to output if there are no threshold crossings 

IF NALARM = 0 THEN GOTO STEP 7 

4. Compute the estimated signature location and signature 
amplitude scale factor for each alai^m. 

FOR I - 1 TO NALARM 

FOR K = L(I) TO N-NH 
IF A(K) = 0 THEN 

KE - K (NEGATIVE-GOING 

THRESHOLD CROSSING) (5.2-9) 


K = N-NH (EXIT LOOP) (5.2-10) 

NEXT K 

TMAX = -1E38 (5.2-11) 

KB = L(I) (5.2-12) 

FOR K = KB TO KE-1 

IF Y(K) > TMAX THEN TMAX = Y(K) (5.2-13) 

LL = K (5.2-14) 


NEXT K 
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L(I) = LL (LOCATION OF I-TH DETECTION) (5.2-15) 

A(I) s Y(LL)/SNR (SCALE FACTOR FOR I-TH DETECTION) 

(5.2-16) 

CR(I) = l/(2n*FM*Y(LL)) (C-R LOWER BOUND) (5.2-17) 

NEXT I 

5. Coapute EN, the expected nuaber of false slaras caused 
by aodeled noise in the residual altiaeter data. 

EN s (N-2«NH)FAR (5.2-18) 

6. Coapute A£, the standard deviations of the errors of 
the estiaated signature aaplitude scale factors 

AE s 1/SNR (5.2-19) 

7. Output: NALARM, (L(k): k » 1,2, . . . ,NALARN} , 

(A(k): k - 1,2,. . . ,NALARH} , EN 

8. End 
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6 . 


GEOSTROPHIC -VELOCITY ESTIMATION ALGORITHM 


6.1 THEORY 

Nearly geostrophic boundary currents, such as the 
Gulf Strean in the western North Atlantic, produce character- 
istic signatures in altimeter data when the satellite subtracks 
intersect the current. Figure 1.3-1 depicts a satellite sub- 
track crossing a current at an angle 6. At position x along 
the subtrack, the dynamic sea-surface height H(x) is related 
to the cross- track component of the geostrophic velocity V^<x) 
by the equation 


V^(x) = - f ^ (6.1-1) 


f = 20sin^ = Coriolis parameter 
0 = earth's rotational velocity 
# ~ north latitude 
g = acceleration of gravity 

For the hyperbolic -tangent current signature H(x) 
described in Section 1.3, the cross- track velocity is 

Vc<x) = \ j * ■ sech^O x/Vg) (6.1-2) 

8 

A ~ amplitude of dynamic height change 
= width of current (90% height change) 

= W^/sin0 = width of signature when sin 0/0 
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6 = track angle with respect to current velocity 
The naxiBUB cross- track geostrophic velocity is 

iVcW • Vc«» = rH; 

Equation 6.1-3 indicates that the BaxiBUB cross-track 
velocity can be coBputed froB estiBates of the aBplitude A and 
width of the boundary-current signature. The required esti- 
Bates of A and are coBputed froB a track of residual altiBeter 
data in two steps. 

In the first step, the residual altiBetry are processed 
with a bank of five Batched filters; each filter is optiBized 
for a different width of signature (V^ * 50 » 60, 75, 100, 150 
kB). The particular filter that produces the largest scaled 
test statistic in response to the boundary current 

is identified; its value for V. is the BaxiBuB- likelihood es- 

tiBate V for the saBple of data being processed. 

8 

A 

In the second step, the BaxiBUB- likelihood estiBate A 
of the signature aBplitude is c<»puted as 

A = A-A, (6.1-4) 

A = signature aBplitude used in designing the 
Batched filter 

A 

A = estiBated aBplitude scale factor based on 

* Batched- filter output 

K * JY(t)l^^SNR 

SNR = Batched- filter ras signal -to-noise ratio 

The estiBated BaximuB cross-track velocity is then 
coBputed as 
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6,2 SPECIFICATIONS FOR THE THRESHOLD DETECTION ALGORITHM 

The geostrophic-velocity estimation algorithm is named 
GVE. The formal specifications for this algorithm are given in 
the following. 

NAME: GVE 

PURPOSE: Estimate the , location and magnitude of the maximum 
cross- track velocity in a geos trophic current. 

REQUIRED EXTERNAL 

ALGORITHMS: ACOVAR, DESIGNMF. CONVOLVE, and DETECTOR 

ACOVAR - autoregressive modeling algorithm 
DESIGNMF = matched- filter design algorithm 
CONVOLVE = matched- filter implementation algorithm 
DETECTOR ~ threshold detection algorithm 
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INPUTS: N, {D(k); k « {LAT(k); k«l»2 N} 


N * nuttb«r of data in the realdual altlsatry tlM 

D(k) = k^^ saaple in the residual altisetry tise aeries 

LAT(k) s north latitude of k^^ saaple in altlaetry 
tine series 


OUTPUTS: KL, V, LE, VE 

NL * estiaated location of aaxiaua geostrophic 
velocity (data saaple nuaber) 

V = estiaated aaxiaua cross-track geostrophic 
velocity (aeters per second) 

LE • Craaer-Rao lower bound on ras error in estiaated 
location NL of the aaxiaua geostrophic velocity 
(saaples) 

VE ~ standard deviation of the error in the velocity 
estiaate V, due to aodeled noise (aeters per 
second) 


ALGORITHM: 

1. Use the algorithas ACOVAR, DESIGNMF. CONVOLVE, and 
DETECTOR to test if a geostrophic current signature 

is in the residual altiaetry data (D(k); k * 1.2.... .M}. 
For this test, input to DESKWMF the tanh aodel current 
signature defined in Section 1.3 with the following 
paraaeter values: A * 0.5 a*. « 75 ka. 

If a geostrophic current is not detected, then set 
VsO. NLsO, L£s0. VE«0 and skip to Step 7. 

2. Truncate the data set (D(k); k * 1.2.....M} to rMK>ve 
the detected geostrophic current signature. Process 
the truncated data set with ACOVAR to coapute an 
awtoregressive (AR) aodel. 

3. Use the AR aodel coaputed in Step 2 with the algoritha 
DESIQfMF to design a bank of five aatched filters, 
each aatched to a tanh current signature having a 
different width. Suggested widths are V * 50. 60. 

75, 100. and 150 ka. The signature aaplftudes are 
all set to A s 0.5 aeter. 
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4. Isplcnent the five matched filters frcMi Step 3 by 
using the algorithm CONVOLVE. Filter the residual 
altimeter data (that part which contains the detected 
current signature) five times , once with each of the 
matched filters. 

5. Compare the scaled test statistics Y(t) coming from 
each of the five matched filters in Step and iden- 
tify that output sequence which achieves the largest 
value in the vicinity of the geostrophic current. 

6. Use the algorithm DETECTOR to compute the geostrophic 
current location NL, the Cramer-Rao lower bound LE on 
the rms location error, the signature amplitude scale 

0k A 

factor A., and the rms error AE in the A. estimate, 
s s 

Compute the estimated maximtm geostrophic velocity V 
3 g A 0.5 

V s =» ; W * signature width corresponding 

2 f W * to largest Y(t) in Step 5 

g = 9.81 m/8^ 

f =7.29 X 10*^ • 2 • sin(LAT(NL)) 

CMBpute the standard deviation VE of the error in the 
velocity estimate V 

VE « V-AE/A^ 

7. Output: NL, V, LE, and VE 

8. End 
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7 . 


SUMMARY 


This report documents the specifications of NOSS algo* 
rithes for ocean current napping and suanariaes the detection 
theory on which the algorithns are baaed. The inputs to the 
algorithns are individual tracks of residual satellite radar 
altineter data frcsi which estinatecl geoid profiles have been 
subtracted. The algorithns are based on the fact that cold* 
core and warn* core current rings and boundary currents can be 
detected by identifying the occurrence of characteristic sea* 
surface height signatures in the residual altineter data. In 
the cnse of nearly geostrophic boundary currents, the cross* 
track c<Msponent of the current velocity can be inferred by 
estinating the alpng* track sea-surface slope fron the altineter 
data and then using the geostrophic equation to conpute the 
velocity. 


Optinal natched filters are used to detect, locate, 
and estinate the anplitudes of generic current signatures in 
the residual altineter data. The algorithns autonatically 
analyze each track of residual altineter data and conpute an 
optinal autoregressive nodel for the noise signal in the data. 
Using this noise nodel, together with a paranetric nodel for 
the deteministic ocean-current signature that is to be detected, 
the algorithn designs a statistically optinal natched- filter 
detector for discrininating between the noise and the signature. 
The detector is optinal in the sense that the probability of 
detecting the ocean-current signature is naxinized for a spec- 
ified probability of false alam (a false alam occurs when the 
randon noise excursions in the altineter data nasquerade as a 
current signature and cause a false detection). The algorithn 
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adjusts the sensitivity of the detector to achieve a specified 
average false-alam rate (e.g., 1 false alarm per 10,000 km of 
data along the satellite subtrack). 

The algorithm for estimating the geostrophic veloci> 
ties of boundary currents eiq>loys a bank of five matched* filter 
detectors: each filter is matched to a different width for the 
current signature. The algorithm determines that signature 
width which is most probable (given the available altimeter 
data) and computes a maximum* likelihood estimate of the current 
signature asq>litude. From this information, the algorithm 
estimates the maxianim along* track slope of the sea surface and 
uses the geostrophic equation to compute the estimated geo* 
strophic velocity. The rms accuracy of the velocity estimate 
is also computed by using the Cramer*Rao lower- bound on the 
variance of the estimated signature amplitude. 

The develofmient of the algorithms is documented in 
Ref. 12. This reference also suwaarizes the results of veri* 
fication tests in which the algorithms are used to detect cold* 
core current rings and to estimate geostrophic velocities with 
SEASAT'A altimeter data and Marsh-Chang geoid estimates (Ref. 13). 
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